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Abstract 

Data from Earth sensor assemblies (ESAs) often are used in the attitude determination (AD) for both spinning and 
Earth-pointing spacecraft. The ESAs on previous such spacecraft for which the ground-based AD operation was 
performed by the Flight Dynamics Division (FDD) used the Earth scanning method. AD on such spacecraft 
requires a model of the shape of the Earth disk as seen from the spacecraft. AD accuracy requirements often are 
too severe to permit Earth oblateness to be ignored when modeling disk shape. Section 2 of this paper rexamines 
and extends the methods for Earth disk shape modeling employed in AD work at FDD for the past decade. A new 
formulation, based on a more convenient Earth flatness parameter, is introduced, and the geometric concepts are 
examined in detail. It is shown that the Earth disk can be approximated as an ellipse in AD computations. 

Algorithms for introducing Earth oblateness into the AD process for spacecraft carrying scanning ESAs have been 
developed at FDD and implemented into the support systems. The Tropical Rainfall Measurement Mission 
(TRMM) will be the first spacecraft with AD operation performed at FDD that uses a different type of ESA — 
namely, a static one — containing four fixed detectors D, (/ = 1 to 4). Section 3 of this paper considers the effect of 
Earth oblateness on AD accuracy for TRMM. This effect ideally will not induce AD errors on TRMM when data 
from all four D, are present. When data from only two or three D/ are available, however, a spherical Earth 
approximation can introduce errors of 0.05 degrees to 0.30 degrees on TRMM. These oblateness-induced errors 
are eliminated by a new algorithm that uses the results of Section 2 to model the Earth disk as an ellipse. 

1. Introduction 

Data from Earth sensor assemblies (ESAs) often are used in altitude determination (AD) for both spinning and Earth- 
pointing spacecraft. Past spacecraft with AD operations performed at the Flight Dynamics Division (FDD) used Earth 
scanning ESAs. The Tropical Rainfall Measurement Mission (TRMM) will be the first such spacecraft to use a different type 
of ESA — namely, a static Earth sensor assembly (SESA), built by Barnes Engineering Company, that comprises four 
detectors Z), (/ = 1 to 4) to sense the Earth horizon at four equally spaced points. 

AD algorithms that use ESA data require a model of the shape of the Earth disk as seen from the spacecraft. Frequently, AD 
accuracy requirements are too severe to permit the simple approximation that the Earth is a sphere, yielding a perfectly 
circular disk. In such cases, a more exact disk shape model, which includes the effect of Earth oblateness, is needed. The 
methodology prescribed by Liu' has provided the basis for Earth oblateness and disk shape modeling at FDD since the late 
1970s. 

Techniques for including Earth oblateness in AD and for correcting attitudes in which it has not been included, on either 
spinning or Earth-pointing spacecraft using scanning ESAs, have been used at FDD for many years^’^'*. No analogous 
techniques were available previously at FDD for spacecraft that carry the Barnes SESA. As discussed later, such spacecraft do 
not require a correction for Earth oblateness in AD if the spacecraft, like TRMM, is nominally oriented along the geodetic 
nadir or horizon bisector nadir (rather than along the geocentric nadir) and if data from all four D, are available. However, 
conditions in which data from only two or three of the four D, are present are not infrequent. While attitude (roll and pitch 
angles) can be computed with data from only two or three adjacent D , , the simple AD algorithms previously used in these 
cases permit Earth oblateness-induced AD errors. 


' This work was supported by the National Aeronautics and Space Administration (NASA)ZGoddard Space Flight Center 
(GSFC), Greenbelt, Maryland, Contract NAS 5-31500. 
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This paper is divided into two subsequent sections. In Section 2, Liu’s Earth disk shape analysis is re-examined and extended, 
and some of the concepts are explored in more detail. The structure of Section 2 largely follows that of Reference 1, and some 
of the material presented in Section 2 has appeared earlier in References 5 and 6. In Section 3, a new AD algorithm is 
prescribed for TRMM when data from only three or two D, are present. 

2. Earth Horizon Modeling 
2.1 The Earth Spheroid (ES) 

Previous AD work at FDD involving nonspherical Earth modeling approximated the Earth as an oblate spheroid prescribed 
as follows; 

x^+/+z^[l-/p=a^ (1) 

where a is the equatorial radius, and /is the flattening. The figure specified by Eq. (1) is called the Earth spheroid (ES) here. 
It was noted recently by Challa^, however, that using a different flattening parameter a related to /as follows 

a = -l+(l-/)' (2) 


is convenient because it facilitates the use of vector methods. Assuming/^ 1/298.257, Eq, (2) yields a = 1/148.379; thus the 
Earth’s a is approximately 2/ 

Eq. (2) enables Eq. (1) to be written as 

-a^ (3) 


where f is a unit vector from the geocenter, O, along the North pole axis, and r is the geocentric radius vector to an 
arbitrary point on the ES. 

We view Eq. (3) as a constraint on r , which we call the ES constraint. Our immediate goal is to define, on the ES, the Earth 
horizon contour as seen from a spacecraft, S. This requires specifying one more constraint on r , and this topic is addressed 
in the following sections (2.2 and 2.3). 

2.2 The Horizon Spheroid (HS) 


Define the following scalar field 


/(r) • z]^ 


and let g{F)be its gradient 


g(r ) -2 [r + a (F • z) f] 


( 4 ) 

( 5 ) 


We view Eq. (4) as specifying an infinite family of oblate spheroids, each with constant equatorial radius sqrt (/(F)) and all 
elements of the family having flattening a and center at O. Let E[r) denote the spheroid at the lip of F and let T(F) be the 
tangent plane to E(f) at F . Then g(F) is orthogonal to r(F) . 


Let 5 be a vector from O to a spacecraft, S. Let h{r) be a vector from S to F . Obviously, 


h{F) = F 


- s 


( 6 ) 


Our second constraint on F is that S, and hence h[F ) , lie in T(F) . Then, 


^(F)»g(F)=0 (7) 

Eq. (7) is the first form of the second constraint on F . Satisfaction of this constraint implies that the lip of F is a horizon 
point of E{F) . A more explicit form of Eq. (7) can be obtained by introducing Eqs. (5) and (6) into it to produce 
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( 8 ) 


r • [f “ ?] + a [F • f ] [f • (r - 5)j = 0 

To transform Eq. (8) into a form that clearly shows the shape of the figure defined by the second constraint, add (J/2) • {sj!) 
to both sides and factor. For simplicity, then introduce v = r - sf2 . The result is 

- - r -,-|2 

V • V + a\v • zY = — • — + a f«— (9) 

^ J 2 2 [ 2 ] 

Eq. (9) shows that the second constraint requires the tip of r to lie on an oblate spheroid, called the horizon spheroid (HS) 
here. Similarly, we call the second constraint the horizon spheroid constraint. The HS is centered at s/2 and, like the ES, 

has flattening a. The equatorial plane of the HS is parallel to that of the ES, and its equatorial radius is the square root of the 
term on the right side of Eq. (9). Points O and S both lie on the HS. 

2.3 The Horizon Plane (HP) 

Eqs. (3) and (8) constitute a pair of constraints on r which, when solved for F , define a closed curv^c, called here the horizon 
contour (HC), which is the ES horizon seen from .S' and specified on the ES. Geometrically, the HC is the intersection of two 
oblate spheroids; namely, the ES and HS, as illustrated in Fig. 1. 


Horizon Spheroid (HS) 



Earth Spheriod (ES) 

Figure 1. Earth Spheroid, Horizon Spheroid, and Horizon Plane 


The geometry and the solution for r , however, can be simplified by replacing the HS constraint by a new constraint formed 
from the ES and HS ones. To develop this new constraint, first note that Eqs. (4) and (5) show that r •g{F) = 2/(F) and 

Eqs. (6) and (7) indicate F • g{F) = s • g{r) . Hence 

?*g(^) = 2/(r) (lO) 

Eq. (10) is merely another form of the HS constraint. Eqs. (3) and (4) show that the ES constraint can be written as merely 

/(F) = «' (11) 

The initial equation for the new constraint is produced by inserting Eq. (11) into Eq. (10) to obtain trivially 
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S • g(r) = 2a^ 


( 12 ) 


For a more useful form of the new constraint, first introduce Eq. (5) into Eq. (12) and manipulate to obtain 

iff? = a^ 


where 


w = ^[5 + of cos A'] 


(13) 

(14) 


In Eq. (14) we have factored s into s =ss and introduced Z' as the co-latitude angle between 5 and f . Using m = mm, 
Eq. (14) indicates that 

m = j^l + 2a(l4-0.5a)cos^ A'j ^ (15) 


Define the following vector n along m : 


where 


Eq. (13) then can be written as 


n=nrh 


n - fm 


P •fh~n 


(16) 

(17) 

(18) 


Eq. (18) is the second form of the new constraint. Analysis of Eq. (18) shows it constrains the tip of r to lie on a plane whose 
normal is along m and whose distance from O is n. This plane and the new constraint are called the horizon plane (HP) and 
HP constraint here. An alternate derivation of Eq. (18) is presented in Ref 5. 

Some later-needed properties of the HP will be delineated before proceeding. Fig. 2 shows the geometry. With the exception 
of ^ , all vectors on this figure lie in a plane normal to the HP. O' is the tip of the perpendicular from O to the HP, and O" 
is the intersection of s with the HP. The position vectors of O' and O” with respect to O are given by h and s^s , 

respectively. The distance between O” and S is given by ^ 2 ; thus, .y, + ^2 ^he term b is the distance between O' and 

O" . The unit vectors E and N are directed along local East and local North, respectively, at O' and define the HP. The 
triad {E,N , m) defines a reference frame centered at O' . The vector m makes angles a and with f and f , 
respectively. Note that is a small angle, since the ES is very nearly spherical; the size of has been exaggerated on 
Fig. 2 for clarity. 


A 

N 



Figure 2. Horizon Plane Geometry and Variables 
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The following important relations can be derived using the material presented above: 


coscr = m*z = cosA' 

(19) 

sincr = [<y/w] sin>^' 

(20) 

C0S6’^ = f •/« = \sjm^ Jl -i-acos^ A'j 

(21) 

sin 6:^ = [,y/w] asinA' cosA' 

(22) 

E - z X s/sinA' 

(23) 

mx E ~ ^5/(wsin/l')j|f^l + acos^ A'^- 5(l + a)cos2' 

(24) 


/? = 


cos/l' sinA’ |l + acos^ A'^ 

= |a^/5'||l + acos^ 2'j 
1 -h acos^ A' - Jl + a cos^ A'^ 


Note that Eqs. (21) and (15) yield 

cos^^ «1 

Thus, the angle between .v and m is given by 


sin^ A' cos^ A’ 


^asinA’ cos A' 


(25) 

(26) 

(27) 

(28) 
(29) 


except for the cases A' =0 and A' = njl . Thus, the geoccnter-to-spacecraft vector ? generally differs from the normal to the 
HP by an angle of the order of the flattening factor a. 

It will be shown in Section 2.4 that the locus of the horizon points in the HP is an ellipse, and it will then be shown that the 
spacecraft views this ellipse slightly obliquely. 

2.4 The Horizon Ellipse (HE) 

The HC now is specified by Eqs. (3) and (18) as the intersection of a spheroid and a plane; namely, the ES and the HP, as 
illustrated in Fig 1, To solve these equations, express r as 

r = n+w (30) 


with vv • « = 0 . Thus, vv lies in the HP, with its base at O' . Eq. (30) thus requires r to satisfy the HP constraint. To further 
restrict r to satisfy the ES constraint, insert Eq. (30) into Eq. (3). After manipulation, we obtain 


w • w + 


a ^•z -I- 2a[/7 • f][vv • f -n^n~a^n•z^ 


z? 


(31) 


Eq. (31), in conjunction with Eq. (30), specifies the HC. 

Eq. (31) can be transformed into a simpler form by specifying the position vector of the horizon point relative to O" rather 
than relative to O' . Letting q be this new position vector and noting that 0"0’ = bN as is evident from Fig. 2, we thus 
have 

w = q- Nb (32) 
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with q , like w , restricted to the HP. 

When Eq. (32) is inserted into Eq. (31), it can be shown, with much work, that a pair of terms cancel out, yielding the 
following simpler result: 

+ -oi^E (33) 

where 

=a[5//w]^ sin^ A' (34) 


•« -a[rt •f]^ + j 

= 1 + acos^ A'-(dj/j)^J|l + crcos^ A'J 


(35) 


Eq, (33) shows that the shape of the HC in the HP is an ellipse, called the horizon ellipse (HE) here. The center of the HE is 
at O” , its flattening is a and its semi-major axis is along E . The HC’s elliptical shape and the important fact that its 
center is on s were noted in Ref 5. 


2.5 The Earth Disk (ED) 

We use the term Earth disk (ED) here to denote the figure of the Earth as seen from S'. The next task is to specify the ED’s 
shape from the preceding material on the HP and HE. Fig. 3 shows the geometry. 


s 



Figure 3. Earth Angular Radius p Geometry 


The previous section (2.4) showed that the HP is not quite orthogonal to s . So we will define a new plane, called here the 
rotated horizon plane (RHP), which is orthogonal to s . The RHP is formed by rotating the HP about E through angle . 

Let g = be a vector on the HE at an azimuth angle from E . Lei p = pp bQ the vector in the RHP corresponding to 
q . Let y/ p be the azimuth angle of p from E . p is formed by projecting q along the vector h onto the RHP, where h is 
the horizon vector from S to the tip of q . Obviously, the locus p{^ ^ in the RHP is not an ellipse centered at O ” . 
Designate the small angle between p and q as (^. The following important equations can be derived: 

p = 9COS<'[l+(9/52)sin4'] ' (36) 

sin4'= (37) 


234 


(38) 

(39) 


C0$4' = Cj 

COS^^ =Cj COS£:^ COS^^ 


E| =[l-sin^ £-mCOS^ ¥p 

q = aE^\+aE&ya^ ¥q\ 


- 0.5 


(40) 

(41) 

(42) 


Let p be the angle at S between h and the unit geocentric nadir vector c = -s . The shape of the ED will be specified by the 
locus of p[¥^ p) around it. It turns out that the simplest result is obtained using the arcctn function: 


p - arcctn (^2 / p) 


(43) 


Insert Eq. (36) into Eq. (43), and use Eqs. (37) through (42) and earlier equations as necessary. After much work, the 
following important result ensues: 


p = arcctn < 1 + a cos^ A.' 


a sin /I' cos /I' sin^^/ -f 


ctn Pc 1 + a[i 


+ k2 sin^ \!/ 


.)i" 


(44a) 


where 


= cos^ X* 

1 + sec^ (l + acos^ 

(44b) 

k^ = sin^ X' |l + acos^ X' sec^ 

(44c) 

> 

Oc =arcsin(a/5) 

(45) 


Pc above is merely the usual angular radius, from S, of a fictitious spherical Earth with radius a and center at O. 

Eqs.(44a) through (44c) were derived earlier in Ref 6 using a different approach. Ref 6 shows, by analytical and numerical 
means, that Eqs. (44a) through (44c) are equivalent to Eq. (4-24) of Ref 1. 

2.6 The Earth Disk Ellipse Approxiniation 

The Earth’s a is small enough to permit Eqs. (44a) through (44c) to be linearized in it with acceptable accuracy. The result 
is 

ctnp = ctnp^ |l + 0.5a |cos^ X' sec^ Pc sin^ Pc + 2sin/l'cos/l' tanp^ sin^^ + sin^ sin^ ^/^j] 

Eqs. (44) and (46) specify disk angular radius p relative to the geocentric nadir . We wish to generalize Eq. (46) to specify 
disk radius relative to lines close to z^. . Fig. 4 shows the geometry on the unit celestial sphere (UCS). Point P here lies on 
the extension of z^ onto the UCS. The reference origin is to be moved from P south through small angle 6' to a new point 
by transforming Eq. (46) from variables p, to new ones y/ For now, e is assumed to be arbitrary, but small 
enough to permit small angle approximations and ruthless linearizations. Fig. 4 then yields 


Ctn p~c\npc\\-^s sec p^ esc sin ^ J 

(47) 

sin^ y/ p- sin^ ^^^1 - 2 6*ctnp^ cos^ y/ ^ esc j 

(48) 

sin ^ p = sin ^ |l - £• ctn p^. cos^ ^ ^ esc ^ j 

(49) 
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Insert Eqs. (48) and (49) into Eq (46) and equate to Eq. (47). This produces an equation in the desired terms e, if/ c- 
However, it is not evident how to solve it for p^ in closed form. 

Horizon Point 



Figure 4, Origin Transiation on Unit Ceiestiai Sphere (UCS) 

The difficulty will be handled by now restricting the work to a single, highly important case. Specifically, € now will be 
limited to the value that causes the two terms involving itself and sin respectively to cancel. The subscript HB, for 
horizon bisector, will be used to denote this ^:and terms pertaining to it. The resulting equation for is 

£jfg = a cos A' sin A' sin^ (50) 

It should be realized that is the angle between 2 ^ and a line, , called here the horizon bisector nadir. An alternate 

name, the apparent nadir, has been used by Flatley^ who performed the original work on the concept. The resulting equation 
for cinpfjQ when Eq. (50) is used is 

ciiip^Q =cln/7^. 1 +0.5a |cos^ 2' sec^ sin^ + sin^ sin^ (51) 

Solving Eq. (51) for pjjg and linearizing in a produces 

P}fB ~ E V'/zb] (52a) 

where 

- Pc - 0.5a cos^ A' sin^ p^ sec p^ (52b) 

= O.S{al Pc) sin Pc cosp^ sin^ X* (52c) 

Eqs. (52a) through (52c) show that the ED can be approximated by an ellipse, centered at p^^ with flattening and semi- 
major axis a^ along local East. This result was derived earlier in Ref. 6 by different means. A MATLAB study, described 

in Ref 6, indicated that for the Earth a, the ellipse approximation is accurate to about le-5 degrees. It should be realized, 
however, that the ellipse exists only on an abstract mathematical plane, not on a physical one. Geometrically, the figure of 
the Earth seen from S can, of course, best be considered to exist on the UCS, not on a plane. 

3. Application to TRMM Attitude Determination 
3.1 Introduction 

Fig. 5 shows the nominal TRMM SESA FOV geometry on a circular ED. The basic outputs of the SESA are the four positive 
penetration angles X, between the Earth IR horizon and the detector D, bases Spacecraft attitude is specified by a 2 x 1 
vector ^of the angular deviation of the SESA and spacecraft body axis 2 ^ = z^ from the HB nadir line z^^ There are two 

main steps in the AD process. First, the are transformed into an SESA frame, S, attitude vector using the 
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algorithms discussed below. Second, is transformed into a spacecraft body, B, frame vector by a rotation 

p =RP . Our cunent interest is only in the first of these two operations. The (j> and 6 are roll and pitch angles, 
determination of yaw is done separately and requires a separate sensor. 



Figure 5. Nominal SESA FOV Geometry on Earth Disk (ED) 

3.2 The Standard SESA AD Algorithms 

The standard algorithms for computing from the X, imply a spherical model of the Earth When all four X, are 
available, we use 

<j>s=Q.5[X^-X^\ (53) 

(95 =0.5[X2 -X,] (54) 

If only three X, are available, either Eq. (53) or Eq. (54) will be unusable, requiring a new algorithm for the affected angle. 
For example, if X 3 is missing, then is computed via 

(!>5=X4-0.5[X, +X 2 ] (55) 

and if X| is missing, 0^ is computed via 

es = X2-Q.5\X^+X^] (56) 

Analogous algorithms are used when X^ or X^ is missing. The function of the ral terms in Eqs. (55) and (56) is to 
compensate for deviation of altitude from that which yields all X, = 0 in the nominal spherical Earth altitude condition. 

When the X, from two adjacent D, , say D, and Dj , are missing, we use 

<i>s=X,-X^y^ (57) 

^S~^2~^AVa 

with 
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~ ^ r\^ ~ ^r\ (59) 

where h is actual altitude, is a reference altitude, is the nominal spherical Earth X^ at h ,^ , and is a scale factor. 

3.3 Proposed AD Algorithms to Compensate for Earth Oblateness 

Section 2.6 demonstrated that to first order in or, the ED is an ellipse with center P^g on Zgg . Therefore, due to symmetry. 
Earth oblateness will not introduce an error into the attitude computed by Eqs. (53) and (54). Numerical studies using the 
software package MATLAB showed that, in fact, this result does not depend on the first order approximation in a. Thus, 
Earth oblateness does not induce an AD error when all four A', are available and used. 

Our MATLAB studies, however, also showed that Earth oblateness induces a significant AD error on TRMM when only 
three or two X, are available and Eqs. (55) through (59) are used. This error varies from 0 to 0.05 degrees in the 3 - X, case 

and from 0.05 to 0.30 in the 2 - X, one. To eliminate the AD errors in the 3 -X, and 2 - X, cases, we have proposed that 

Eqs. (55) through (59) be replaced by a new algorithm based on the elliptical ED model from Section 2.6. The main steps of 
this algorithm are as follows. 

(1) Compute the azimuth angles y/ of the D, relative to local East via 

Wd,b (60) 

where t^ris the spacecraft yaw angle, is the constant azimuth angle of each D^ relative to xg , and y/ ^ is the spacecraft 
velocity angle relative to local East, y/ ^ can be computed using 

y/^ = njl - arctan (cos //cos (eu + /)) (61) 

where /. m, and/ are orbit inclination, argument of perigee, and true anomaly, respectively. 

(2) Use the y/^ to compute the ED radii /?£,_ at the D, using Eq. (52). Then compute the nominal penetration angles 
^iNOM via X^oj^ = ¥d^ where 4 is the detector mounting angle on Fig. 5. 

(3) Compute the and/or 6 g that cannot be computed by Eq. (53) and/or Eq. (54) via 

~ ^iNOU (62) 

~ ^ j ~ (63) 


where X, and X ^ are the proper available penetration angles. Theory indicates this algorithm should almost perfectly 
eliminate Earth oblateness-induced AD errors in the 2 - X, and 3 - X; cases. 
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